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In general neutron stars in binaries are spinning. Recently, a new quasi-equilibrium approximation 
that includes a rotational velocity piece for each star has been proposed to describe binary neutron 
stars with arbitrary rotation states in quasi-circular orbits. We have implemented this approximation 
numerically for the first time, to generate initial data for neutron star binaries with spin. If we choose 
the rotational velocity piece such that it equals the Newtonian rigid rotation law, we obtain stars 
with fluid 4- velocities that have expansion and shear of approximately zero, as one would expect for 
quasi-equilibrium configurations. We also use the new approach to construct and study initial data 
sequences for irrotational, corotating and fixed rotation binaries. 
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I. INTRODUCTION 

Binary neutron stars are at the intersection of two of 
the most fascinating topics in astrophysics: Gamma ray 
bursts and gravitational wave astronomy. Binary neutron 
star mergers (together with black hole neutron star merg- 
ers) have been proposed as potential engines for short 
duration gamma ray bursts [lH8|. These are likely gen- 
erated in the massive accretion disks around the merger 
remnant: a larger neutron star or a black hole. In ad- 
dition, binary neutron star systems are one of the most 
promising sources for gravitational wave detectors such 
as LIGO d, OH, Virgo p], H3 or GEO [3. Several 
of these detectors have been operating over the last few 
years, while several others are in the planning or con- 
struction phase [lij ]. During the inspiral regime, when 
the two stars are still well separated they can be well 
approximated by post-Newtonian theory. Later, when 
the stars get close, their matter distributions eventually 
merge together to form a single differentially rotating ob- 
ject. Depending on the total mass, the two progenitors' 
spins, the equation of state and the strength of mag- 
netic fields this object can either promptly collapse to a 
black hole, or form a hypermassive neutron star. The hy- 
permassive neutron star is supported against collapse by 
differential rotation. It can survive for many dynamical 
timescales, while angular momentum is gradually trans- 
ported from the inner to the outer parts. Eventually the 
hypermassive neutron star will also collapse to a black 
hole surrounded by a massive torus, that is more mas- 
sive than in the prompt collapse case. Such systems 
could supply the energy required for a short gamma ray 
burst IML% 

In order to make predictions about the last few orbits 
and the merger of such systems, fully non-linear numer- 
ical simulations of the Einstein Equations are required. 
To start such simulations we need initial data that de- 
scribe the binary a few orbits before merger. The emis- 
sion of gravitational waves tends to circularize the or- 
bits [la Jl6l]. Thus, during the inspiral, we expect the 
two neutron stars to be in quasi-circular orbits around 



each other with a radius that shrinks on a timescale 
much larger than the orbital timescale. This means 
that the initial data should have an approximate heli- 
cal Killing vector In general these neutron stars will 
be spinning. In the case of the double pulsar PSR J0737- 
3039 [ljj the spin period of the faster spinning star will be 
Pa(£) = 27ms at merger |l8[ and should thus not be ne- 
glected. In the quasi-circular regime the orbital timescale 
will be much shorter than the spin precession time scale, 
thus we can assume that the spins are approximately 
constant. 

To incorporate these ideas and to construct numerical 
initial data for binary neutron stars with arbitrary spins 
and masses we will use an approach introduced in [l8| . 
In this approach the stars are given spin by choosing a 
rotational velocity for each star. In this way it is possible 
to construct stars with both rigid or differential rotation. 
In equilibrium we of course expect the stars to be rigidly 
rotating such that the expansion and shear of the fluid 
4-velocities of each star vanish. We find that this can 
be achieved (to good approximation) by setting the rota- 
tional velocity of each star equal to the Newtonian rigid 
rotation law. 

We also construct and discuss initial data sequences 
for irrotational, corotating and fixed rotation binaries. 

Spin will have a noticeable effect on the inspiral and 
merger of the binary if the spin period is within an or- 
der of magnitude of the orbital period. Since the orbital 
period is on the order of a few milliseconds during the 
late inspiral, we expect interesting spin effects for spin 
periods on the order of a few dozen milliseconds or less. 
To date we have observed only ten binary neutron star 
systems, thus it is not clear yet how likely such spin pe- 
riods are. Some population synthesis models [l9[ suggest 
that radio observable pulsars in neutron star binaries do 
have a distribution of spin periods that extends down to 
about 15ms. Other population synthesis models for pul- 
sars [20| come to similar findings. However, since such 
models involve many parameters that are used to de- 
scribe sometimes poorly understood physical processes it 
may be too early to draw definite conclusions about what 
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spin periods can be expected in neutron star binaries. 
One parameter that may be of particular importance and 
that illustrates this uncertainty is the magnetic field de- 
cay timescale Td- In the models it is assumed that the 
magnetic field of each star decays exponentially on this 
timescale. Since neutron stars spin down due to magnetic 
dipole radiation, magnetic field decay can have important 
effects for the spin down rate and thus the expected spin 
periods before merger. Unfortunately the value of Td is 
controversial. In the first model mentioned [l9[ the mag- 
netic field decay timescale has to be r d ~ 5 My in order to 
fit observations, while in the other model [2fJ one needs 
r d ~ 2000My. 

Throughout we will use units where G = c = 1. Latin 
indices such as i run from 1 to 3 and denote spatial in- 
dices, while Greek indices such as p run from to 3 and 
denote spacetime indices. The paper is organized as fol- 
lows. Sec. [IT] lists the General Relativistic equations that 
govern binary spinning neutron stars described by perfect 
fluids. In Sec. IIIII we briefly describe what algorithm we 
use to numerically implement these equations. We then 
present some numerical results in Sec. IIVI In particular 
we describe how one should choose the rotational velocity 
of each star. We also present particular examples in the 
form of several initial data sequences. We conclude with 
a discussion of our method in Sec. [Vj Some derivations 
are relegated to the appendices O and [5] 



II. BINARY NEUTRON STARS WITH 
ARBITRARY ROTATION STATES 

In this section we briefly describe the equations gov- 
erning binary neutron stars in arbitrary rotation states in 
General Relativity. These equations were derived in [l8j . 

We use the Arnowitt-Deser-Misner (ADM) decomposi- 
tion of Einstein's equations and describe the gravitational 
fields (i.e. the 4-metric g^v) in terms of the 3- metric 7^, 
lapse a, shift /3 l and the extrinsic curvature K lJ . Wc 
further assume that the neutron star matter is a perfect 
fluid with stress-energy tensor 



[p (l + e) + P]«V +Pg> MV . 



(1) 



Here po is the mass density (which is proportional the 
number density of baryons), P is the pressure, e is the 
internal energy density divided by po and w M is the 4- 
velocity of the fluid. Assuming a polytropic equation of 
state 



D 1 + 1 A 
P = k Po 



(2) 



and defining the specific enthalpy 

h = 1 + e + P/p (3) 
Po, P and e can all be expressed in terms of h. We find 



where we have used the abbreviation 
q = (h-l)/(n+l) 



(5) 



The fluid 4-velocity u M is expressed in terms of the 
3- velocity 



hHu" 



(6) 



which in turn is split into a irrotational piece D l <j) and a 
rotational piece w l 



(7) 



where Di is the derivative operator compatible with the 
3-metric 7y. 

In order to simplify the problem and to obtain ellip- 
tic equations we make several assumptions. The first is 
the existence of an approximate helical Killing vector £ M , 
such that 



£t9iiv ~ 0. 



(8) 



We also assume similar equations for scalar matter quan- 
tities such as h. For a spinning star, however, is 
non-zero. Instead we assume that 



0, 



(9) 



so that the time derivative of the irrotational piece of the 
fluid velocity vanishes in corotating coordinates. We also 
assume that 



7f £v*w„k 0, 



and 







(10) 



(11) 



which describe the fact that the rotational piece of the 
fluid velocity is constant along the world line of the star 
center. 

These approximations together with the further as- 
sumptions of maximal slicing 



Hi 



K tJ = 



and conformal flatness 



7y = lp 4 Sij 
yield the following coupled equations: 

V' 5 



D 2 %\) - 



32a ; 



■(LB)^{LB) ij +2tt^p = 0, 



(12) 
(13) 

(14) 



Dj(LBy j - (LB) 13 Dj ln(aV>~ b ) - I6iraip 4 f = 0, (15) 



Po = K- n q n 



P 

e 



K' n q n+1 
nq, 



(4) 



D 2 (aip) — aip 



7jA 
32o7 



{LB) t3 {LB) l3 +2^ i {p + 2S) 



= 0, 
(16) 
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Di 



(£>V + O-p mx u (/? l +e) =0, (17) 



where we have introduced 



and 



h = ^L 2 - (Drf + Wi){D*4> + w l ). 



(18) 



Here (LB) lJ = D l B J + D^> B % - p^D k B k , D t = d u and 
we have introduced 



(19) 



p = a 2 [po(l + e) + P]u u° - P, 
f = a[po(l + e) + P}u°u (u l /u + (3 l ), 
S ij = [ Po (l + e) + P]u u (u i /u + l3 i )(u j /u + P j ) 
+Pl l \ (20) 



o _ Vh 2 + (D^ + w t )(D^ + w*) 
an 

2 _ b+^/b 2 - 4a 4 [(A</> + w 4 )«/] 2 
L ~ 2a^ ' 

6 = [{e+p i )D i cj ) -Gf + 2a 2 {D i 4> + w i )w\{21) 

where we sum over repeated spatial indices irrespective 
of whether they are up or down and where C is a constant 
of integration that in general has a different value inside 
each star. Below we will denote these two values by C+ 
and C- . Notice that in an inertial frame the approximate 
helical Killing vector has the components 



(l,-fi[ 



Mix 1 



(22) 



Here x l CM denotes the center of mass position of the sys- 
tem, and n is the orbital angular velocity, which we have 
chosen to lie along the £ 3 -direction. 

The elliptic equations ([T3]) . (fi~5j) . (fTB]) and (HH) above 
have to be solved subject to the boundary conditions 

lim ip = l, lim B l = 0, lim aip = 1 (23) 

r— >oo r— »oo r—tao 

at spatial infinity, and 

(-DV)Apo + w l D iPo = hu°(p l + C)D iP a (24) 

at each star surface. Note that the rotational piece of the 
fluid velocity w % can be freely chosen. 

Once the equations ([14]), ([15]), ([16]), ([17]) and ([18]) arc 
solved we know h (and thus the matter distribution) and 
the fluid 3- velocity ^u 1 via Eq. ([7]). The 3-metric is ob- 
tained from Eq. (|13[) and the extrinsic curvature is given 
by 



1 



2?/> 4 a 



(25) 



Notice that Eq. (|T5|) can also be written as 



In h 



hu°) 



r = u\ 



i3 i + e + 



hu° J a 2 hu° (ahu )' 



.(27) 



lnr + ln(-C), 



(26) 



Our approach, while more general, has certain simi- 
larities with the approach in [2l], [22[, hereafter BS. For 
example our Eq. ([7]) reduces to Eq. (52) of BS if we as- 
sume w % = rj^ 1 . Furthermore, using the approximations 
in Eqs. ([9]), ([TO]) and (fTT|) . we are able to find the first 
integral (fT8|) of the Euler equation. This means that 
within our approximations the Euler equation has van- 
ishing curl. That is precisely what has to hold if the 
approach in BS (which converts the Euler equation into 
an elliptic equation by taking its divergence) is to suc- 
ceed. The problem with the BS approach is that it simply 
assumes £^u^ — instead of our Eqs. ([9]), ([TO]) and ([TT]) . 
This leads to extra terms in the Euler equation, which 
cause a non- vanishing curl. Thus, the BS approach is not 
entirely consistent since it requires zero curl of the Euler 
equation, while at the same time it uses an Euler equa- 
tion with non-vanishing curl. Furthermore, the boundary 
condition given by BS for their new elliptic equation has 
to be imposed at the star surface. Since the location of 
the star surface is another unknown function one needs 
an equation or an algorithm to determine this location. 
No such algorithm is given by BS. The usual iterative 
approach where we simply search for the surface where 
h = 1 in each iteration (see e.g. Sec. IHIj) will not work, 
because the boundary condition given by BS ensures that 
h = 1 occurs at the surface where we impose their bound- 
ary condition. So in this kind of iterative procedure the 
star surface would always remain at the location of the 
initial guess for the surface. 



III. NUMERICAL METHOD 

To construct initial data we have to solve the elliptic 
equations ([11]), (JTU), ^ and (TTT]) together with the al- 
gebraic equation (TT8l) . This set of equations has a similar 
structure as for the well known case of irrotational neu- 
tron star binaries, which has been solved before [23M3TI ] . 
In this work we will use the SGRID code (32l[33 |. which 
uses pseudospectral methods to accurately compute spa- 
tial derivatives. We use the same decomposition into 6 
domains as was used in [33[ for the case of corotating 
neutron star binaries. In this approach one star center is 
put on the positive and the other on the negative x-axis. 
Then complicated coordinate transformations are used 
to transform from Cartesian like coordinates (x, y, z) to 
new coordinates (A, B, <p). Here the coordinates A and B 
both range from to 1, and ip is a polar angle measured 
around the x-axis. The actual coordinate transforma- 
tions are different in each domain. This results in two 
domains that cover the outside of each star (including 
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spatial infinity) for either x > or x < 0. These two do- 
mains touch at x = 0. Since they contain spatial infinity 
it is trivial to impose the boundary conditions in Eq. ([23)) . 
The coordinate transformations contain freely specifiable 
functions cr±(B, ip) so that one can always make the inner 
domain boundaries coincide with the star surfaces. The 
inside of each star is covered by two more domains. One 
of these stretches from the star surface up to a certain 
depth inside the star. The other covers the remainder of 
the star interior. The elliptic equations (fT4"j) . (fT5|) and 
(|16p need to be solved in all domains, while the matter 
equations (fTTf and (IT51) are solved only inside each star. 
Two of the domain boundaries always coincide with the 
neutron star surfaces so that it is straightforward to im- 
pose the boundary condition (|24|) for <fi at each star sur- 
face. Notice, however, that Eq. (fTT)) and its boundary 
condition in Eq. (|24|) do not uniquely specify a solution 
<j>. If 4> solves both Eqs. (fTTJ) and ([M]) (j> + const will be 
a solution as well. In order to obtain a unique solution 
we modify the boundary condition by adding the volume 
integral J star 4>dV over the star interior to the left hand 
side of Eq. (j2~4")l . Furthermore, in the domains covered by 
the A, B, ip coordinates we impose the following regular- 
ity conditions along the x-axis: 



o, d s * + d a d v d v v = 0, 



(28) 



where stands for either ip, B l , a or <f), and s = 
\f y 2 + z 2 is the distance from the x-axis. 

In order to solve the elliptic equations (fT4"]l . (|T5|) . (fll))) 
and ([17]) we need a fixed domain decomposition. How- 
ever, the location of the star surfaces (where h = 1) is 
not a priorily known, but rather determined by Eq. (TT51) . 
For this reason we use the following iterative procedure: 

1. We first come up with an initial guess for h in 
each star, in practice we simply choose Tolman- 
Oppenheimer-Volkoff solutions (see e.g. Chap. 23 
in [35|) for each. For the irrotational velocity po- 



tential we choose <fi = fX^c* 



and , 



b CMJ 



where x\ 



C* 



■CM 



is the center of the star and the center of 



mass. We choose the initial orbital angular velocity 
according to post-Newtonian theory. 

2. We then evaluate the residuals [i.e. the i 2 -norm 
of the left hand sides of Eqs. dTtj, flU]), (HO) and 
(fT7|)] . If the combined residual is below a prescribed 
tolerance we are done and exit the iteration at this 
point. 

3. If the residual of Eq. (H) is larger than 10% of 
the combined residuals of Eqs. (|T4"|) , (p~5|) and (fl"6|). 
we solve Eq. (fT7)) for cj>. We then reset <fi to </) = 
0.2(j) e ii + 0.8(f) o id, where 4> e u is the just obtained 
solution of Eq. (|17l) and tfi id is the previous value 
of <f>. 

4. Next we solve the 5 coupled elliptic equations (fl4l . 
(1T51) and CLU) for * eH = (ip,B\a) eH . We then set 
* = (i/j, B\ a) to = 0.4*di + 0.6* oid . 



5. In order to also solve Eq. (|18|) we need to know the 
values of the constants C± in each star as well as 51 



and Xq M . We first determine the star centers x\ 



c*± 



by finding the maximum of the current h along the 
x-axis. Since the location of each star center is 
given by <9i/i| x i_ = 0, Eq. (|21)|) [which is equivalent 

to Eq. (0] yields 



di In 



/3 ,; + C + 



hu° J 



•6 + 



hu°) 



= -25ilnr (29) 



Note that f3 l + ^ is a function of f2 and x\, M . 
We now update SI and x x CM by solving Eq. 
for f2 and x x CM so that the star centers Xc*± re- 
main in the same location, when we update h ac- 
cording to Eq. (fT5|) or Eq. (|2"6| . For this reason 
Eq. (|29p is sometimes referred to as force balance 
equation. One noteworthy caveat is that we eval- 
uate the derivative of lnT in Eq. ([29)1 for the SI 



and Xq M before the update. Since we iterate over 
these steps this approximation does not introduce 
any errors. However, it is essential for the overall 
stability of our iterative procedure. 



6. Next, we use Eq. (fT8| to update h in each star, 
while at the same time adjusting C± such that 
the rest mass of each star remains constant. This 
update is numerically expensive because the do- 
main boundaries need to be adjusted (by chang- 
ing <j±(B,<p)) such that they remain at the star 
surfaces, which change whenever h is updated. 
When we adjust a±(B,ip) it can be helpful for 
the stability of the overall iteration to filter out 
high frequency modes in a±(B,(p) and to impose 
8b^±{B, <p)|b=o,i = 0. The latter keeps the stars 
from drifting away from the x-axis during the iter- 
ations. 

7. Finally we go back to step 2. 



IV. NUMERICAL RESULTS 

We have implemented the method described above in 
the SGRID code [32H34| . In this section we present some 
numerical results using this code. All our results are 
presented in units where G = c = k = 1 and we only use 
n = 1 in the polytropic equation of state. 



A. Choice of rotational velocity piece w* 

As already mentioned we hold the rest masses fixed 
during the iterations described in Sec. IIIII Similarly we 
need to fix the rotational velocity piece w l that gives 
rise to the spin in each star. Ideally we would like to 
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choose w l such that the expansion and shear of the fluid 
are approximately zero. As we show in Appendices [X] 
and [Bj this is true if the velocity in the corotating frame 
= u»/u° - ^ is of the form 



v y 



yi = e ijku] j \x h 



£.(*)]• 



(30) 



where cc^ is the location of the star center, which could 
be defined as the point with the highest rest mass density 
Po or as the center of mass of the star. Thus our task 
is to choose w % such that Eq. (f30f holds. In |l8[ we 
have speculated that a good choice might be DiW 1 = 0, 
which can be rewritten in terms of the derivative operator 
Di = di as DiW 1 = if we introduce 



w 



(31) 



It is clear that 



w 



f(\x" 



U C* I 



ijk 



U C*J> 



satisfies DiW 1 = for any function f(\x n 



(32) 
|) that 

only depends on the conformal distance from the star's 
center. We have tested the simplest case of / = 1, and 
find that the V 1 that results from this choice after solving 
all equations is similar to Eq. (|30[) but with a position 
dependent u> J , so that we have differential rotation in 
the star which leads to non-zero shear. Another simple 
choice is 



(33) 



Numerically, we find that this choice results in a V 1 that 
is very close to the form in Eq. (|3T)|) . In Fig. [T]we show 
results for an unequal mass system where both stars have 
w l as in Eq. J3p) with u/ = (0,0,0.025). Since V i varies 
linearly with x l we see that it is of the form in Eq. (I3TJ1) . 
The slope of V % corresponds to Cj % » (0,0,-0.015), so 
that we obtain io % w ft 1 + uj 1 , which would hold exactly 
in Newtonian theory. 



B. Initial data sequences 

In order to test our method we have performed sim- 
ulations with different rotation states. In Figs. [2] and 
[3] we show how the ADM mass Madm and the ADM 
angular momentum Jadm vary as a function of the or- 
bital angular velocity f2 for a binaries with rest masses 
moi = 0.1461 and TO02 = 0.1299. Note that increasing fl 
corresponds to decreasing separation. As we can see our 
numerical results (pluses, stars and crosses) approach the 
expected post-Newtonian results (taken from [36 39] ) for 
non-spinning point particles (solid line) for small Q. In 
fact we find that the results for irrotational (w l = 0) 
stars (marked by stars in Figs. [5] and [3]) agree very well 
with post-Newtonian non-spinning point particle results 
for all ft we have investigated. On the other hand, for 
corotating binaries (pluses) we obtain larger Madm and 
Jadm values, especially for larger f2, which is expected 




FIG. 1: The Cartesian x- and y-components of V 1 in the 
orbital plane of a the star with rest mass moi = 0.1461. 
The other star has rest mass mo2 = 0.1299. The separa- 
tion between both star centers is D — 5.2885, which cor- 
responds to an orbital angular velocity of Q = 0.03928. 
Both stars' rotational velocity piece is given by Eq. ()33|) with 
u/ = (0,0,0.025). We see that both V x and V v vary linearly 
and thus obey Eq. ((30)) . 



M, 



post-Newtonian 
NR: corotating 
NR: irroational 
x-x NR: co = 0.025 




0.2568 



0.2566 



0.07 



FIG. 2: The ADM mass for a binary with rest masses moi = 
0.1461 and mo2 = 0.1299 as a function of orbital angular 
velocity. Shown are results for post-Newtonian point particles 
(solid line), and three different numerical results (NR) for 
corotating stars (pluses), irrotational stars (marked by stars), 
and a case where both stars have spin (crosses) with u 1 = 
(0,0,0.025). 



because to maintain corotation the stars have to spin 
more for higher Q. Finally, we have also investigated 
the case of a binary where both stars have the same 



constant rotational velocity w' 



uj 3 (x k 



b C* 



) with 



uj 1 = (0,0,0.025). This uj 1 corresponds to a spin pe- 
riod of 14ms [for k = 0.018m 5 /(kg s 2 )]. In Figs, d] and 
[3] this case is denoted by crosses. Since here the stars 
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NR: corotating 
*-* NR: irroational 
x-x NR: co = 0.025 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 

FIG. 3: The ADM angular momentum for the same binaries 
as in Fig. [2] 



always have rotational velocity we obtain larger Madm 
and J adm values than in the irrotational case for all il. 
However, since the rotational velocity is always the same 
we get less Madm and J adm than in the corotating case 
(pluses) for large f2, and more Madm and J adm than 
in the corotating case (pluses) for small Vl. 



V. DISCUSSION 

Realistic neutron stars in binaries are spinning. From 
observations of millisecond pulsars we know that these 
spins can be substantial enough to influence the late in- 
spiral and merger dynamics of the binary. The spins 
might also influence the lifetime of an angular momen- 
tum supported hypermassive neutron star that can form 
after merger. With more angular momentum we expect 
the hypermassive star to survive for longer. This can 
have important consequences for both the gravitational 
waves emitted by the system and also for the likelihood 
of a gamma ray burst. 

The purpose of this paper is to numerically implement 
and test a new method for the computation of binary neu- 
tron star initial data with arbitrary rotation states. This 
method is derived from the standard matter equations 
of perfect fluids together with certain quasi-equilibrium 
assumptions. We assume that there is an approximate 
helical Killing vector £ M and that Lie derivatives of the 
metric variables with respect to vanish. We also as- 
sume that scalar matter variables such as h or po have 
Lie derivatives that vanish with respect to £ M . However, 
since the Lie derivative of the fluid velocity is non- 
zero for generic spins, we split the fluid velocity u^ L into 
an irrotational piece (derived from a potential <p) and a 
rotational piece w % , and assume that only the irrotational 
piece has a vanishing Lie derivative (see Eq. ((HP) with re- 
spect to Furthermore we know that the spin of each 



star remains approximately constant since the viscosity 
of the stars is insufficient for tidal coupling [4(| . To in- 
corporate this fact, we use Eqs. (fTUj) and (fTTj) which are 
based on the assumption that w l is constant along the 
star's motion described by the irrotational velocity piece 

From these assumptions we obtain the elliptic equa- 
tions (flip, (flip. (fT5P and (HZP together with the algebraic 
equation (fTSp . The specific enthalpy h in Eq. (TT51) deter- 
mines the shape of the star surfaces. Since our elliptic 
solvers work on a fixed domain decomposition where the 
star surfaces have to coincide with domain boundaries we 
solve this mixture of elliptic and algebraic equations by 
iteration. In each iteration we first solve the elliptic equa- 
tions for a given h and then use the algebraic Eq. (jT5J) 
to update h. The stability of this iterative procedure is 
improved if we do the following. (A) We typically do 
not take the ip, B l and a and 4> coming from solving 
Eqs. (flip. ([I5p. (fTBT) and (Q7P as our new fields. Rather, 
we take the average of this solution and the values from 
the previous iteration step as our new fields. In this way 
ip, B l and a and <fi change less from one iteration step 
to the next. (B) We use the force balance condition as 
given in Eq. (|29| to update f2 and x x CM . 

For each iteration we also need to specify rotational 
piece w % of the fluid velocity. We have found that the 
choice in Eq. (l33l) works very well in the sense that after 
numerically solving all equations it leads to a velocity in 
the corotating frame of the form V = Co x r as in Eq. (|30|) . 
In Appendices [5] and [B] we show that this form of V 
results in a fluid 4- velocity with an expansion and shear 
that are approximately zero, as we would expect for stars 
in equilibrium. This means that we have found a simple 
way to generate initial data for neutron star binaries that 
ensure that the star are spinning and without differential 
rotation. 

We also compare initial data sequences for irrotational, 
corotating and fixed rotation binaries. We find that our 
method yields reasonable results. All sequences approach 
post-Newtonian results for large separations. And the 
binary sequences with fixed rotation yield higher Madm 
and J adm than corotating configurations for large sep- 
arations, and lower Madm and J adm than corotating 
configurations for small separations. 
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Appendix A: Expansion, shear and rotation of the 
fluid 



If m m denotes the 4- velocity of the fluid we can split its 
derivatives into 



VvUi 1 = -9P liU + cr liV + u) ia ,-a l j,Uv, (Al) 



where 



(A2) 



and where the expansion, shear, rotation and acceleration 
are defined as 



(A3) 



= P» P v v V^tV) - -6P^, (A4) 



^ = P»P V V V kV ] 



and 



If the 4- velocity is of the form 

it" = /u", 



(A5) 
(A6) 
(A7) 



where / is any scalar function, it immediately follows 
(from P liv u v = 0) that 



9 = fP^V M u v , 



(A8) 



<v = fPg K v (M ^, - -ejv, (A9) 



and 



(A10) 



For our purposes it is often convenient to write the 
4- velocity as 



(AH) 



where in an inertial frame the helical Killing vector is 
given by Eq. (|22[) and u° = —n^u^ 1 /a. Note that £° = 1 
leads to V° = 0. From 



dt 71° ? 



(A12) 



we see that V* can be interpreted as the velocity in the 
corotating frame. 

The fact that Vr^^) = together with Eqs (|A8|) and 
HI yields 



9 = u u Pi""V M Vu, 



(A13) 



and 



= uOpH'pfv^V^ - -QP^. (A14) 



Thus we see that if we have V ' (pV v ) ~ 0, the expansion 
and shear are approximately zero. 

Appendix B: Expansion and shear if V = U) X r 



Let us assume that V 1 = (0, V 1 ) is given by Eq. ([SO]). 
Let us further assume that V 1 is small in the sense that 



V i ~ O(e) 



(Bl) 



with t «1. We know (e.g. from post-Newtonian theory) 
that the shift is also small. For simplicity we assume that 

P l ~ 0(e). (B2) 

Using V( AI V r t/ ) = £vg^ v j2 and Eq. (fT3)) we then find 



V (0 U ) - -V l d ia + 0(e 2 ) 
V(o^) = 0(e 2 ) 



(B3) 



Recall that the Newtonian potential near a mass mi is 
given by 



m i m 2 I -, , x — x Cif {t) 
D 



(B4) 



where is a second mass at a distance D in the x- 
direction. Since tp ocU and a oc U we find that 



V 1 ^ ~ y'fta ~ 0(e 3 ), 



(B5) 



where we have used the virial theorem and set ^ 
0(e 2 ). 

From u^u^ = — 1 we obtain 



tt = ± +O(c 2) 

a 



(B6) 



which leads us to 

P^ = 7^ + 0(e 2 ). (B7) 
From Eqs (|aT3)) , ([ATI)) . ([Bljl. (|B5]l and JB7| we see that 



6 = C7 MV = 0(e 3 ) 



(B8) 



if y is of the form uj x r as in Eq. ([3H)l . This means that 
expansion and shear for this particular V % are smaller by 
a factor of 0(e 2 ) than for a generic V 1 . 
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